#include<../data_tools/usefull_routines.c>

int main(int argc, char *argv[])
{
  int i,j,nbines,ntypes,nfiles,liminf,limup,account;
  double *radio,rmin,rmax,deltar,rini,rfin;
  double vcprom;
  size_t *ind;
  char outfiles[1000];
  particulas *diskpart;
  
  printf("\nComputing curve of rotation for:\n");
  nbines = atoi(argv[1]);
  nfiles = atoi(argv[2]);
  int npart[nfiles];
  char *infiles[nfiles];
  ntypes = 0;
  for(i=3; i<nfiles+3; i++)
    {
      infiles[i-3] = argv[i];
      npart[i-3] = contadorParticulas(infiles[i-3]);
      ntypes = ntypes+npart[i-3];
      printf("%d particles from file %d\n",npart[i-3],i-2);
    }
  printf("%d particles\n",ntypes);
 

  diskpart = (particulas *)malloc((size_t)ntypes*sizeof(particulas));
  if(diskpart == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  radio = (double *)malloc((size_t)ntypes*sizeof(double));
  if(radio == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  ind = (size_t *)malloc((size_t)ntypes*sizeof(size_t));
  if(ind == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  liminf = 0;
  printf("Reading circular velocities\n\n");
  for(i=0; i<nfiles; i++)
    {
      FILE *typefile;
      typefile = fopen(infiles[i],"r");
      printf("%s\n",infiles[i]);
      limup = liminf+npart[i];
      for(j=liminf; j<limup; j++)
	{
	  fscanf(typefile,"%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
		 &diskpart[j].id,
		 &diskpart[j].pos[0],&diskpart[j].pos[1],&diskpart[j].pos[2],
		 &diskpart[j].vel[0],&diskpart[j].vel[1],&diskpart[j].vel[2],
		 &diskpart[j].mass,
		 &diskpart[j].pot,
		 &diskpart[j].posc[0],&diskpart[j].posc[1],&diskpart[j].posc[2],
		 &diskpart[j].velc[0],&diskpart[j].velc[1],&diskpart[j].velc[2],
		 &diskpart[j].mang[0],&diskpart[j].mang[1],&diskpart[j].mang[2],
		 &diskpart[j].mr,
		 &diskpart[j].r,
		 &diskpart[j].vc,
		 &diskpart[j].epsilon);
	  radio[j] = diskpart[j].posc[0];
	  ind[j] = j;
	  // printf("%d\n",j);
	  //if(radio[j]<0.00001) printf("si\n");
	 
	}
      liminf = limup;
      fclose(typefile);
    }
    
  gsl_sort_index(ind,radio,1,ntypes);
  
  printf("Done.\n");
  printf("\n");

  printf("Printing curves of rotation\n");

  rmin = log10(radio[ind[0]]);
  //rmax = log10(radio[ind[ntypes-1]]);
  rmax = log10(20.0);

  deltar = (rmax-rmin)/nbines;
  rini = rmin;

  FILE *vcfile;
  if(nfiles!=1)
    {
      vcfile = fopen("total_curve_rotation.dat","w");
    }
  else
    {
      sprintf(outfiles,"%s_curve_rotation.dat",infiles[0]);
      vcfile = fopen(outfiles,"w");
    }
  for(i=0; i<nbines; i++)
    {
      rfin = rini+deltar;
      vcprom = 0.0;
      account = 0;
      for(j=0; j<ntypes; j++)
	{
	  if((radio[ind[j]]>=pow(10,rini)) && (radio[ind[j]]<pow(10,rfin)))
	    {
	      account = account+1;
	      vcprom = vcprom+diskpart[ind[j]].vc;
	    }
	  // if(radio[ind[i]]<0.00001) printf("si\n");
	  //if(i==npart[0]) fprintf(vcfile,"FILE 1\n");
	  //if(i==npart[1]+npart[0]) fprintf(vcfile,"FILE 2\n");
	  
	}
      vcprom = vcprom/account;
      fprintf(vcfile,"%e %e\n",pow(10,rini),vcprom);
      rini = rfin;
    }
  printf("Done.\n");
  printf("\n");

  free(diskpart);
  free(radio);
  free(ind);
  fclose(vcfile);

  return 0;
}
